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ABSTRACT 



We introduce a method for obtaining the X-ray luminosity function (XLF) 
and the binary-period distribution of populations of high-mass X-ray binaries 
(HMXBs) in the stellar fields {i.e., outside globular clusters) of normal galax- 
ies. We start from standard distributions of the parameters of those primordial 
binaries which are the progenitors of HMXBs, and follow the transformation of 
these ditributions with the aid of a Jacobian formalism as the former evolve into 
the latter through the processes of the first mass transfer and the supernova 
(SN) that follows. We discuss the distributions of the post-SN binaries and the 
HMXBs. We show that our calculated model XLF has a differential slope ~ —1.6 
with a flattening at low luminosities, in excellent agreement with observations. 
The calculated binary-period distribution, which basically has a slightly slop- 
ing plateau-like character at intermediate periods, with a rise to this plateau at 
shorter periods and fall-off from it at longer periods, is in agreement with the 
observed distribution within observational uncertainties. We discuss the physical 
origin of these distributions. We demonstrate that, while the effects of both (a) 
the distribution of the properties of the massive companion in the HMXBs, and 
(b) the primordial orbital distribution and the SN dynamics are important, the 
former appear to be dominant in determining the XLF, and the latter in deter- 
mining the HMXB binary-period distribution. We discuss the possible roles of 
stellar-mass black holes and ultra-luminous X-ray sources (ULX) in the observed 
"universal" XLF of HMXBs. 

Subject headings: binaries: close - stars: neutron - stars: massive - supernovae: 
general - X-rays: binaries - X-rays: galaxies 
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Introduction 



Four decades have passed now since tlie Uhuru discovery of t 



le X-ray binary 



(henceforth XRB) Cen X-3 was first reported (iGiacconi et aL 197ll ). What began in 
the early- to mid-1970s as detailed individual studies of the small number of XRBs 
known at the time, as also pioneering works on the construction of general mod els of 



such binaries ( Lamb et al. 1973 



Prinde k Rees 1972 



Davidson fc Ostriker 19731) has 



now blossomed into a mature, rich subject (for overviews, see, e.g., Shapiro & Teukolsky 
1983; Ghosh 2007), with a large number of XRBs discovered and catalogued today in the 
Milky Way, the Magellanic Clouds, and nearby external galaxies, and the fine details of 
the behavior of each class of XRBs well-recorded from timing and spectral studies with 
the aid of several generations of X-ray observatories of ever-increasing sensitivity and 
resolution. It is thus becoming possible now to do studies of the collective properties of 
XRB popul ations, obtaining statistically meaningful d i stributions of some of their essenti al 



properties fiKim fc Fabbiano 2004 



Grimm et al. 2002 



Grimm et al 2003l: iGilfanov 2004 



Gilfanov et al. 2004a : 



Gilfanov et al. 2004b 



Kim fc Fabbiano 2010 



Liu et al. 2005 



Liu et al. 20061 ) . and exploring the essential physics underlying these distributions, which 
is rooted in the basic scenarios for the formation, evolution, and dynamics of XRBs (van 
den Heuvel 1983, 1991, 1992, 2001), which have been constructed gradually over these four 
decades, and which are widely accepted now. By collective properties, we mean here the 
distributions of X-ray luminosities and binary periods of these XRBs, the distribution of 
pulse periods of those binaries which exhibit periodic X-ray pulses, and perhaps also the 
distribution of their spectral properties at a future date. 

In this series of papers, we explore the theoretical underpinnings of the observed 
distribution of the collective properties of XRBs in normal galaxies. We focus here on 
XRBs outside the globular clusters of these galaxies {i.e., in the stellar field of the galaxy). 
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where the effects of encounters between XRBs and neighboring stars, as well as those 
between XRBs themselves, are generally thought to be negligible, due to the relatively 
low stellar density there. This situation is at the opposite limit to that obtaining in 
dense cores of globular clusters, where such encounters are thought to be dominant in 
determining the formation, evolution, and destruction of XRBs, and so their collective 
properties. We have studied this latter limit of XRBs in globular clusters in an earlier series 
of papers (Banerjee & Ghosh 2006, 2007, 2008). The current study may thus be regarded 
as complementary to the earher one. We shall confine ourselves here to the distributions of 
X-ray luminosity and binary orbital period, since the most detailed distributions available 
are on these properties, and also since these properties are most readily associated with 
binary evolutionary characteristics. To understand the distributions of properties like pulse 
periods of X-ray pulsars and spectral parameters of XRBs, one must also take into account 
the details of the accretion torque on the neutron star, and of the X-ray emission processes. 
We defer these projects to the future, noting that with ~ 140 X-ray pulsars already known, 
the pulse-period distribution may be amenable to such a future study. 

The crucial simplification that makes a straightforward approach feasible in the limit 
that we study here is that the absence of any significant influence on a given XRB of either 
the stellar background or other XRBs imphes that each XRB can be thought of as evolving 
in isolation, following the well-studied and now widely-accepted scenarios for the formation 
and evolution of individual XRBs, starting from primordial binaries generated in normal 
star-formation activity in a galaxy (see, e.g., van den Heuvel 2001 and the references 
therein). This, in turn, implies that we need only have a knowledge of the distribution of 
the essential parameters of these primordial binaries, and evolve these distributions through 
the essential processes that occur during the evolution of a primordial binary into a XRB. 
The task becomes even simpler when only initial and final states matter in the above 
processes, since we need only keep track of these states to carry out the transformation 
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of the relevant probability distribution, and so connect the primordial distribution to the 
XRB distribution in a one-to-one correspondence. This is, in fact, the case for massive or 
high- mass X-ray binaries (henceforth HMXBs). On the other hand, when following one or 
more of these evolutionary processes requires the knowledge of the entire evolution during a 
particular process, and not just of the initial and final states, the calculation becomes more 
laborious (although still straightforward in principle). This is the case for low-mass X-ray 
binaries (henceforth LMXBs). 

In this first paper of the series, we focus our attention on HMXBs containing neutron 
stars. A particular circumstance that helps in this case is that, because of the short 
evolutionary timescales of the massive companions of the neutron stars in HMXBs, the 
timescale for the entire evolution from the primordi al binary stage to the HMXB stage (as 



also the operational lifetime of the HMXB) is short ( Ivan den Heuvel 200ll ) compared to the 



timescale on whic h the star-formation rate (which determine s the creation rate of primordial 



binaries) evolves (IMadau et al. 1998 



Ghosh &: White 200ll ). As a consequence, the entire 



process we study here can be viewed as happening at a constant star-formation rate, so that 
the star-formation history of the galaxy adds no extra complication. Figuratively speaking, 
it is as if we are taking a "snapshot" of the above processes confined entirely to one epoch 
(the current one, z = 0, for the HMXB populations of the Milky Way and other local 
galaxies, and an appropriately earlier epoch for the HMXB populations of spiral galaxies at 
significant redshifts). 

This argument does not, of co urse, apply to LMXBs , since their evolutionary timescales 



are comparable to or longer than (Ivan den Heuvel 19911 ) those on which the star-formation 



rate evolves. Accordingly, we have to explicitly take into account the evolution in the 
star-formation rate when we apply our scheme to LMXBs, as we shall describe in later 
papers in the series. Thus, for example, for describing the LMXB populations of local 
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galaxies, we shall need to use not only the star-formation rate at 2; = but also the 
rates at z > 0, i.e., the star-formation history of the galaxy. This difference between the 
evolutionary characteristics of HMXB and LMXB populations is, of course, closely related 
to the earlier statement about the initial-to-final state correspondence for HMXBs versus 
the additional role of intermediate states for LMXBs. 

In the following sections, we first describe the distributions of the essential parameters 
of the primordial binary populations that we adopt for our study of HMXBs, following 
well-established scenarios and norms in the literature. In particular, we adopt the standard 
initial mass function (henceforth IMF) and the standard log-uniform distribution (also 
known as Opik's law) for the orbital separation, as is done in the literature. We then 
describe how t hese parameters chang e in the essential processes that occur as a primordial 



binary evolves (Ivan den Heuvel 20011 ). The first change occurs when the primary evolves, 
fills its Roche lobe, and transfers its envelope to the secondary, leaving behind its He core. 
The second change occurs when this He core finishes its further evolution, and explodes in 
the supernova (SN) that creates the neutron star. We give explicit relations connecting 
the binary parameters in the initial and final states of these two processes. Using these 
relations, we show how the standard rules of probability transformation enable us to 
derive the distributions of the post-SN binary parameters from the original primordial 
binary parameters through the Jacobian formalism. Next, we consider what happens when 
the massive companion to the neutron star evolves off the main sequence, becoming a 
giant /supergiant, and driving a strong stel lar wind, from which the neutron star accretes 



matter, thus turning on the HMXB phase fivan den Heuvel 20011 ). 



We adopt standard stellar-wind models from the literature and show how we can derive 
the distributions of both the X-ray luminosities Lx and the orbital periods Port of the HMXBs 
from the distributions we obtained above. We then compare our calculated distributions with 
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and Porb (ILiu et al. 2005 



the obser ved distributions of Ly (IGrimm et al. 20021 : iGilfanov 2004 : iGilfanov et al. 2004al ) 



Liu et al. 20061 ) reported in the literature. We show that the main 



feature of the observed Lx-distribution, viz., th e power-law be havior with a differential 
slope dN/dLx ~ —1.6 over a wide range of Lx (IGilfanov 20041 ) is reproduced well in our 



calculated distribution over the luminosity range covered by neutron-star HMXBs. At 
the lowest luminosities, below the range covered in the above observational work, our 
calculated distribution shows a shallower rise in dN/dLx with decreasing Lx- We show 
that this feature is consistent with recent observations (Shtykovskiy & Gilfanov 2005a,b) 
of the Magellanic Clouds. For the Porb-distribution, we show that the general features of 
the observed distribution are in reasonable agreement with the results of our calculations 
within observational uncertainties. 

Subsequent sections of the paper are arranged as follows. In SecJH we give a brief 
overview of the formation of neutron-star HMXBs from primordial binaries and their 
subsequent evolution, and we describe the standard distribution of primordial binary 
parameters that we adopt. In SecJHl we describe how the binary parameters evolve 
through the processes of (a) the first mass transfer in the system and (b) the subsequent 
supernova of the He-core of the primary. In Sec J4.1l we introduce the Jacobian formalism 
for calculating how the distribution of binary parameters transforms as these parameters 
evolve as described in the previous section. In SecJH we discuss the distribution of the 
essential parameters of the post-SN binary. In Sec|5|, we detail the evolution of the post-SN 
binary into a HMXB, summarizing the essentials of the stellar and stellar-wind models 
that provide these details, and indicating the further transformation of distributions needed 
at this point. In SecJHl we discuss our calculated HMXB distributions of luminosity {i.e., 
the X-ray luminosity function or XLF) and binary period, compare them with observed 
distributions, and discuss the effect on our model calculations of varying essential inputs 
like the initial mass function (IMF), stellar wind models, and so on. Finally, in SecJTJ we 
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discuss our results from various angles, including issues about primordial binary parameters, 
roles of Be-star binaries and black holes, physical origins of some essential features of the 
HMXB distributions, and then summarize our conclusions and future outlook. Details of 
some of our calculational methods are summarized in Appendix A. 



2. HMXB Formation and Evolution 

Scenarios for the formation and evolution of HMXBs have been studied in detail since 
the early 1970s, and a "standard" picture has now emerged, which is very well-document ed 



(Ivan den Heuvel 1983 



van den Heuvel 1991 



van den Heuvel 1992 



van den Heuvel 200ll ) 



thus making it unnecessary to recount it here in detail. Very briefly, one starts with a 
primordial binary of two massive stars. The more massive of the two {i.e., the primary) 
evolves faster, ends its main sequence life and expands rapidly to become a giant, the 
whole process occurring on a timescale of ~ 10^ — 10^ years. The giant overflows its Roche 
lobe and starts transferring its hydrogen envelope to the secondary, at which point there 
are two possibilities, depending on the ratio of the thermal timescales of the two stars, 
which, in turn, depends on their mass ratio. The first possibilty arises if the mass ratio is 
not so extreme (quantitative arguments are given in the next section) that the above two 
timescales are not different by more than an order of magnitude, say. The secondary can 
then accept the entire mass transferred by the primary. This leaves behind only the He-core 
of the primary, and the secondary becomes so massive after assimilating the H-envelope of 
the primary that it in fact turns into the more massive component of the system. This is 
the channel through which HMXBs form. The second posibility arises if the mass ratio is 
so extreme that the thermal timescales are disparate by more than the above amount, in 
which case the secondary is unable to accept the mass transferred by the primary, and this 
mass forms a common enevelope (CE) surrounding the secondary and the He-core of the 
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primary. This is the standard channel through which LMXBs are thought to form, which 
we do not discuss further in this paper. 

In the HMXB channel, the system now consists of a He-core plus a massive star, 
which was the secondary earlier, but is now the more massive component. The He-core 
finishes the rest of its evolution rapidly (in ~ 10^ — 10^ years), and explodes in a supernova, 
leaving behind a neutron star which will become the main X-ray emitter in the HMXB. 
The massive companion to the neutron star now finishes its main-sequence evolution in 
~ 10^ years, expands to become a supergiant, and drives a strong stellar wind. Accretion 
from this wind by the neutron star generates X-rays, and the system thus turns on as a 
HMXB. This HMXB phase lasts ~ 10^ — 10^ years, after which the massive companion 
fills its Roche lobe. The subsequent mass-transfer is very large and rapid, which leads to 
considerable mass loss from the system, formation of a common envelope (CE) surrounding 
the neutron star and the He-core of the evolved massive companion, and a total extinction 
of the X-ray source. The outcome of this CE evolution, viz., the formation of a compact 
binary consisting of the neutron star and the He-core of its former massive companion, 
which, after the SN of this He-star, ultimately leads to the formation of a double-degenerate 
system consist ing either of an eccent ric double-neutron-star binary or of two runaway single 



neutron stars (Ivan den Heuvel 200 ll ). lies outside the scope of this work. 



In recent years, variations, modifications and additions to the above standard scenario 
have been considered. For example, it has been suggested that CE evolution of sufficiently 
wide systems during the first mass transfer may still lead to tigh t post-CE systems which 



would contribute to HMXB production (see (ILinden et al. 20091 ) and references therein) 



and that these systems might remain tight even after the SN explosion, turning promptly 
into HMXBs with unevolved companions. Further, under certain ot her circumstances, C E 



evolution might lead to tightly-bound HMXBs with He-rich donors (ILinden et al. 20091 ) 



-lo- 



in this first look at the problem, we shall not consider these fascinating possibilities, but 
rather confine ourselves to the standard picture, since our aim here is to assess the viability 
of our approach in the simplest testing ground, before attempting further refinements. 



2.1. Primordial Binary Distribution 

We consider now the distribution of the essential properties of the primordial binaries 
from which HMXBs evolve according to the standard scenario sketched above. The 
progenitor primordial binary system is described by three essential properties, namely, the 
mass of primary (Mp), the mass of the secondary (M^), and the orbital separation (flo)- An 
equivalent description is in terms of Mp, the mass ratio q = Ms/Mp, and ao- We use the 
second description throughout this work, following the custom of recent literature on the 
subject, and briefly mention the correspondence between the two descriptions in Sec J7.1[ 
Also, as is often done in this subject, we take the primordial binary orbits to be circular for 
the purposes of the problem we study here. 

The constraints imposed by the HMXB problem under consideration here limits the 
allowed range of the above parameters. If we restrict ourselves to only neutron-star HMXBs, 
as we do in this paper, then the p rimary mass Mp is restricted to be between roughly 



9Mq and 30Mq (iHeger et al. 20031 ). This means, of course, that those HMXBs which have 
black holes (either stellar- mass or intermediate-mass) as their compact. X-ray emitting 
components are outside the scope of this work. We return to this point in Sec J7.3[ Next, the 
secondary mass Mg is restricted from below by the requirement mentioned above, i.e., that 
its thermal timescale must not be larger than that of the primary by more than one order 
of magnitude, which implies a rough lower bound Ms > 0.3Mp, since the thermal timescale 
goes roughly as the inverse square of the mass in this mass-range. The upper bound on 
Ms is of course Mg < Mp by definition. This means, therefore, that the allowed range of 
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the primordial-binary mass ratio is roughly 0.3 < q < 1 for HMXB formation. Finally, the 
initial orbital separation is restricted from above by the condition of Roche lobe overflow 
when the primary evolves off the main sequence, without which there would be no mass 
transfer, and the stipulated HMXB formation scenario will not operate. It is customary to 
take this upper limit as oq < lO^i?0, which is not very restrictive. 

Consider now the distributions of Mp, q, and Oq we use in this work, in keeping with 
the common practice in the literature. The Mp-distribution is described in terms of a 
suitable Initial Mass Function (IMF), i.e., the probability density JmIM) of a stellar 
mass being in the range M to M + dM. The I MF is widely taken to be of a power-law 
form fM^M) oc M~", the classic Salpeter IM F (ISalpeter 1955|) correspo nding to a = 2.35, 
and a more recent suggestion being a = 2.7 (IKroupa fc Weidner 20031 ). We adopt the 
Salpeter IMF for our main calculations and study the effect of varying the IMF on the 
Lx-distribution in Sec J6.4.2l 

The g-distribution is generally taken as a power-law in q, given by f{q) oc q^, where the 



exponent (3 may depend on the type of XRB pop ulation being studied. 



bllowing the usage i n 



Linden et al. 20091 ) . 



recent works on population studies of HMXBs (iBelczynski et al. 2008 
we adopt a uniform g-distribution, /3 = 0, in our HMXB studies in this paper. (However, 
we do note that our formalism has the provision for handling non-uniform distributions of 
q, which will indeed be used in our LMXB studies described in subsequent papers in this 
series.) 

The distribution of the orbital separation is almost universally taken in the literature 
to be a loguniform on e (i.e., an e qual number of systems in equal intervals of logao), also 



known as Opik's law (jOpik 19241 ). This implies a probabilty density /(oq) of the separation 



being in the range ao to ao -l- dao as /(ao) oc l/ag. We adopt Opik's law in our work. 



In order to construct the total probibility distribution function of primordial binaries. 
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we need only note that the above pieces of the probabihty density function (PDF) are 
independent of each other, subject to the restrictions in the allowed ranges of values detailed 
above. Thus, subject to these restrictions, the total PDF is given by: 

1 M-° 

fprimo{Mp,q,ao) = — (1) 

i V CLq 

Here is a normalization parameter. Note that, although this PDF is mathematically 
defined over much larger ranges of parameters, we are interested only in the allowed range 
for HMXB formation, as explained above. Hence, we set fprimo = over the forbidden 
range, and so choose the normalization parameter that the integral of the PDF over the 
allowed range is unity. Consequently, although this PDF has no explicit dependence on q 
due to the assumed uniform g-distribution for this HMXB study, the allowed range of q, 
i.e., 0.3 — 1, does have an effect on the value of A^. 



3. Evolutionary Changes in Binary Parmeters 

In this section, we summarize the changes in the binary parameters as the system 
evolves through the first mass-transfer phase, and the subsequent SN of the He-core of 
the original primary. As mentioned in SecJH these changes only require keeping track of 
the relations between the initial and final states in each of the above two processes, which 
makes the ensuing transformations relatively simple. 



3.1. First mass transfer 

In keeping with previous work on the subject, we assume that there is negligible mass 
loss from the system during this process, i.e., the mass transfer is conservative, so that 
the H-envelope of the primary is entirely transferred to the secondary, leaving behind its 
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He-core of mass Mp^c given by: 

Mp,, = MoM'/, (2) 

where typical values of Mq and ^ used in the literature are O.O73M0and 0.704 respectively, 
for a metallicity of z ~ 0.03. This is a commonly-used analytic approximation to the results 
of numerical stellar-evolution calculations (Ghosh 2007 and references therein). 

Because the mass transfer is conservative, the binary parmeters Mp, Ms, and a at the 
end of it are related to the primordial binary parameters as: 



Mp = Mp,e, 

Ms = Mp + Ms- Mp^^ , 



(3) 



and 



a — ao 



MpMs 
MpMs, 



(4) 



3.2. Supernova 

The post-mass-transfer binary is detached. The He-core described above evolves 
further and explodes as a supernova (SN), which leaves behind a neutron star of a typical 
mass of 1.4 Mq and blows off the rest of the He-star. This sudden mass loss from the 
system alters the orbital separation and makes the orbit eccentric. We assume in this 
work that the neutron star has the above mass in all cases. Neglecting any effects of 
the expanding supernova ejecta on the secondary, e.g., ablation, the secondary remains 
completely unchanged in this process. Since we assume that the mass of the neutron star 
is always fixed to 1.4 Mq, it can no longer be used formally as a free parameter in the 
description of the post-SN system. Therefore, we choose the orbital eccentricity e as the 
third parameter, so that the post-SN system is described by the companion mass (Mc), the 
semimajor axis (a) and the eccentricity (e). 
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Natal supernova kicks introduce additional changes in the post-SN parameters. It is 
generally believed today that the SN explosion need not have exact spherical symmetry. 
Many suggestions have been given for the physical origin of this asymmetry, e.g., density 
inhomogeneity in the pre-collapse core, anisotropic neutrino emission, unequal momentum 



Scheck et al. 20061 ) 



fluxes in the jet and anti-jet directions, and so on ( iScheck et al. 2004 
It has also been pointed out that a seed anisotropy, once introduced by these mechanisms 
will be enhanced during the hydrodynamic evolution of the explosion and may impart 
a large kick to the neutron star. Proper motion measurement of radio pulsars is used 
as a diagnostic of these natal kicks. Such studies have shown that the distribution of 
the pulsar velocities can be described by a 3D is otropic Maxwellian , in which the kick 
speed V is distribut ed as f{v) oc v'^exp{—v'^/2a'^) (IHobbs et al. 20031 ) . with a = 265 km/s 
( iHobbs et al. 20031 ) . and the direction of the kick velocity is distributed isotropically. 
Earlier works had suggested that the distribution might actually be bimodal, represented 
by two Maxwellians, so th at the overall distribution would be the weighted sum of these 



( lArzoumanian et al. 20021 ). 



Recent works have argued that electron-capture supernovae (ECSN) are likely to 



produce much lower kicks with cr 50 km/s or lower (ILinden et al. 20091 ). while the larger 
value given in the previous paragraph is appropriate for iron-core-collapse supernovae 
(ICCSN). We note here that studies of the proper motions of pulsars involve single pulsars, 
some of which would have come from a single star collapse. For the rest, which were 
obtained by disruption of binary systems, there would clearly be an observational bias 
towards higher kick velocities, since it is these systems which became preferentially unbound 
to yield the single pulsars whose proper motions are studied observationally. However, such 
effects are difficult to account for at this stage of understanding of the problem, and it is 
customary in the literature to apply the inferred distribution of kick velocities directly to 
the X-ray binary progenitors, as we have done here. 
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Effects of isotropic Maxwellian SN-kicks in progenitors of X-ray binaries were studied 
in a pioneering work by Kalogera (1996). For our study liere, we have devised a metliod 
related to tliat described in the above work, but designed specifically for our purposes here 
The method works as follows. If we define our co-ordinate system such that, just before 
the SN explosion, the line joining the two stars is along the x-axis and the neutron star 
is moving in the positive z-direction with a keplerian orbital velocity Vorb, then a kick of 
magnitude Vk in the direction n{6, (p) would produce a change in orbital parameters given 



Here the angles 6 and (p are defined in the usual way, M* denotes the total mass of 
the system, and the subscripts i and / denote the initial and final values {i.e., pre-SN and 



Due to the Maxwell-distributed random SN-kicks given to the neutron star, the 
one-to-one correspondence between the pre-SN and post-SN binary parameteres, which 
would have existed in the absence of such kicks, is broken now. Instead, we need to calculate 
a suitable average effect of SN-kicks on the transformation between the pre-SN and post-SN 
parameters, which takes into account the underlying distribution of the kicks. To this end, 
we average eqsj5]over the distribution of v^, upon which the linear terms in VkCosO vanish 
due to isotropy, and the distribution- averaged transformation equations become: 



by: 




post-SN values). All quantities here are in solar units, with v in the units of a/ GMq / Rq . 



— + = mU 

ai \ai afj 
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Here, vl is f ^ averaged over the kick-distribution. For a Maxwellian distribution, it is 
given by: 

/ 2\ ^ j v^exp{-v'^/2(T'^)dv 
J v^exp{—v^/2a^)dv 

Those properties of f| for a MaxweUian which we need for our work are detailed 
in Appendix A. We assume here a combination of two Maxwelhans with values of a 
appropriate for ICCSN and ECSN, as described above. The way for determining the 
appropriate proportions of these two components is described below in Sec J4.2[ For each of 
the two components, the parameter transformations are also given in Appendix A. 



4. Post-SN probability density function 

4.1. Jacobian formalism 

Probability theory provides us with a method of transforming the PDF of a set of 
variable to the PDF of another set of variables, if we know the relations between the these 
two sets of variables. Let us assume that a system is described by a set of variables denoted 
by a vector X. A given specific process transforms this set of parameters to another set of 
(the same number of) parameters denoted by Y. Both the forward transformation Y{X) 
and the inverse transformation X{Y) are mathematically defined. The initial PDF as a 
function of X, fx{X) then transforms as follows: 

= fxiX{Y)) 



dX 



dY 
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Here 



ax 



is called the Jacobian determinant of the inverse transformation. The above 



' dY 



theorem can now be applied to our problem in order to transform the PDF of the primordial 
binaries to the post-SN PDF. 



In order to apply the above formalism to our problem, we need to invert the 
above parameter transformation relations, which is straightforward. The Jacobian of 
this transformation {Jmtr) can also be readily calculated. In the two regions described 
in Appendix A, these transformations can be obtained and the Jacobian (Js„) can be 
calculated explicitly, as detailed in Appendix A. 

Working in either region, the post-SN PDF is given by: 



The PDF of the parameters of the post-SN binary, Mc, e, a, no longer factorizes into 
individual PDFs for these three parameters, unlike the situation for the PDF of primordial 
binaries. This is as expected, since the parameter transformation laws given above mix the 
parameters. An immediate consequence of this is that the PDF for any one of the above 
post-SN binary parameters needs to be obtained by integrating the PDF of Eqnl9]over the 
other two parameters. 

As described above, we use two values of the dispersion in the kick distribution, a, 
appropriate for ICCSN and ECSN. We have also described above the method of calculating 
the distribution in each case. The actual distribution of the parameters with both types of 
SN present is the weighted sum of these two PDFs for the two values of a, and this weighting 
comes from the following considerations. For the range of Mp which leads to ECSN events, 



4.2. Inverse transformation and post-SN PDF 




sn 



(9) 
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we have adopted here the value advocated in recent works flPodsiadlowski et aL 2004 ). 
which is 8-11 Mq. The relative contribution of the ECSN and ICCSN is then given 
by Nec/Nicc = N{Mmin,Mtran)/N{Mtran,Mmax), which depends upon the IMF. Here, 
Mmin = 9M0, Mmax = SOMq define the allowed range of primary masses given earlier, and 
Mtran is the primary mass for transition from ECSN to ICCSN. For a Salpeter IMF, this 
ratio is 0.42 for Mtran = HMq and 0.67 for Mtran = 12Mq. We take this ratio as 0.42 
wherever it occurs in our calculations here. 



4.3. Properties of the post-SN PDF 

In order to appreciate the nature of the post-SN distribution, we calculate and plot the 
individual PDFs for each of these parameters, obtained by integrating fpsn over the other 
two. We now discuss each of these distributions. 



4-3.1. Companion mass 

The distribution of the companion mass shows a slow early rise and then broken 
power-law behaviour with the break occuring at ~ 30Mq, as shown in Fig{T] The power-law 
slope is ~ —2.4, in the mass range 20 — SOMq, whereas above the break point the 
Mc-distribution becomes steeper, with a slope of ~ —5.2. The kink seen at the breakpoint 
of 3OM0 becomes a fold in the bivariate distribution as a function of Mc and a (See Fig. H]). 
This kink appears at the mass which corresponds to the upper limit of the primary mass 
in the primordial binary. The allowed range of the primary mass is 9Mq < Mp < SOM©: 
this constraint was imposed to ensure the formation of neutron star as an end product of 
the stellar evolution of the primary. 

It can be easily seen from the relation given above between the mass of the primary 
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and that of its He-core that the mass range obtained for Mc if we take Mp = 30Mq, 
corresponding to the entire allowed range of Mg, is ~ 30Mq — SIMq. Thus for < SOMq, 
a large range of Mp can contribute, whereas for > 30Mq the allowed range of Mp is 
rapidly cut off as Mc increases, which causes a rapid reduction in the number of systems 
possible in this region. Thus, the effect of the allowed phase-space region is reflected in the 
transition in the power-law index of the Mc-distribution. In the region where large parts of 
the phase space of Mp are allowed to contribute, the slope of the PDF is close to the input 
IMF slope of —2.35. For masses above the kink, only small parts can contribute, resulting 
in a steeper fall. For similar reasons of forbidden phase space, a downturn is expected at 
very low values of Mc, in the mass range ~ (10 — 15) M©. 

Upon inclusion of SN-kicks, the above mass distribution is essentially unchanged in 
the mid-Mc range, except that the kink at ~ 30Mq becomes less prominent. But the 
distribution is cut off more sharply at the lowest and highest masses, the former being 
immediately understandable since kicks tend to preferentially unbind systems with lower 
binding energy. 

The bivariate Mc — a distribution shown in FigH] reveal the complexity introduced 
by the kicks. Although the surface shows similarities to its no-kick counterpart in general 
features, several additional features appear due to Maxwellian-averaged kicks corresponding 
to the two different values of a given above for the two types of SN, namely, ECSN and 
ICCSN. 

4-3.2. Semimajor axis 

Figj2] shows the PDF of log{a), the post-SN semimajor axis on a logarithmic scale. The 
PDF without considering the effect of SN-kicks is flat [i.e., loguniform in a) in the mid-a 
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range, with smooth rise and fall at low and high values of a respectively. This is clearly a 
direct consequence of the assumed flat, Opik's-law distribution of primordial binaries, which 
remains unchanged at the intermediate values of a. The PDF with the SN-kicks included 
clearly shows two additional features. First, there is a general shift of the distribution 
towards wider orbits, entirely as expected, since SN-kicks generally tend to widen orbits. 
Second, the PDF is not exactly loguniform, i.e., flat in the mid-a range now: a slow fall 
with increasing a is observed in this range, before the PDF falls off rapidly at large a, as 
before. This second feature is a consequence of the fact that SN-kicks steadily reduce the 
probability of survival of wider binaries even in the mid-a range, before destroying them 
altogether at large a-values. 

4.3.3. Eccentricity 

In the absence of kicks, the eccentricity introduced by the SN explosion is equal to the 
fractional mass loss form the system. As this fraction is small (~ 10% or less) in case of 
HMXBs, only small eccentricities are introduced by the SN in the no-kick scenario. Figure 
[31 which shows the post-SN eccentricity PDF, confirms this. However, introduction of the 
SN-kicks changes this result completely, making the post-SN systems much more eccentric, 
as expected. Further, the results for ICCSN and ECSN are quite different in terms of 
post-SN eccentricity, also as expected because of the large difference between a in the two 
cases (see above). Due to the small kicks imparted in case of ECSN with independent 
of the details of the pre-SN system (see Sec 13. 2 1 and Appendix A), the eccentricities are 
generally smaller than in the ICCSN case (although larger than in the no- kick case), and 
the e-distribution is spread over a large range. By contrast, for ICCSN the kicks are much 
larger, leading to large eccentricities, as shown in FigJS] Further, in this case vl is almost a 
fixed fraction of the disruption velocity of the pre-SN system for Maxwellian (and possibly 
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also for other similar) distributions, as explained in Appendix A, so that the e-distribution 
is narrowly peaked around a value determined by that fraction. It is likely that more 
detailed calculations would widen this peak somewhat, but no qualitative changes are 
expected. The composite PDF of eccentricity for ECSN and ICCSN naturally shows a 
double-peaked structure, as in FigJHl 

Of course, the e-distribution calculated above refers to the immediate post-SN systems, 
which are not observable as XRBs. These systems subsequently undergo tidal circularization 
rapidly, so that most of them become circular or nearly so by the time mass transfer begins 
and the system turns on as a HMXB. Since eccentricity is an irrelevant parameter at the 
HMXB stage, we integrate fpsn over the eccentricity to obtain only the bivariate PDF as a 
function of and a, which is shown in FigjH 



Post-SN systems evolve over the main-sequence lives of their companions as detached 
systems on a short timescale (~ 10^ years), during which the orbit is tidally circularized. 
Therefore, a bivariate post-SN PDF adequately describes the HMXB system, being given 



Observed collective properties of HMXBs are usually given as distributions of their 
luminosities (L) {i.e., the X-ray luminosity function or XLF) and orbital periods (Pb), 
instead of and a, which we have worked with upto this point. In order to compare our 
results with observations, a further transformation (M^, a) — )■ (L, Pf,) is therefore required 
now. Pf, can be calculated for given values of and a using Kepler's third law. If all the 



5. Post-SN binary to HMXB 



by: 




(10) 
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masses and distances are expressed in the units of the solar values, then Ph in hours is given 
by: 

= 7.72 (11) 

Here MjV5 = IAMq is the mass of the neutron star. To calculate the luminosity as 
a function of Mc and a, prescriptions of the mass-loss rate from the companion and the 
capture mechanism by the neutron star are required. 

The companion expands into a giant /supergiant after completing its main-sequence 
life, and loses mass by driving a strong stellar wind at a rate Mw The neutron star 
captures a fraction of this lost mass. The accretion rate onto the neutron star is then 
given by M = M^x (capture fraction). In this work, we assume that HMXBs are entirely 
wind-fed giant/ supergiant systems, thus ignoring those systems in which a main-sequence 
companion close to its Roche lobe may drive a (weak) atmospheric Roche-lobe overflow. 
This is certainly justified at higher X-ray luminosities, since systems fed only by the latter 
mechanism would be found at the lowest end of the XLF. In a similar vein, we also do not 
include HMXBs with Be-star companions (see Sec J7.2p . since the transient nature of these 
systems with low duty cycles implies a low time-averaged luminosity, which would put them, 
again, at the low end of a time-averaged XLF, which is generally difficult to determine. 

The X-ray luminosity L is then given by the standard stellar-wind model, according to 
which the massive companion drives a wind mass- loss at a rate M^, and with a terminal 

2 

(the capture fraction introduced above) is 
captured and accreted by the neutron star, and the gravitational energy release from this 
accretion generates the X-rays. In order to connect L to the stellar and orbital parameters, 
we first note that L is the rate of release of the gravitational energy of the accreted matter 



wind speed of Vw Of this, a fraction ~ ^^y^ 
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and hence is directly proportional to accretion rate Mas 

^ GMnsM 



R 



NS 



(12) 



which we can express numerically as L^q ^ 1.17M_io, where L^q is L in units of 10^^ ergs/s 
and M_iois Min units of 10""'^'^ Mq/jt. Next, we note that the accretion rate Mis related 
to the wind mass-loss rate M^, and the above capture fraction as 

2 



M = M„ X 



GM. 



NS 



Via 



(13) 



To proceed further, we need models of stellar winds, which we consider next. 



5.1. Stellar wind models 



Models of stellar wind relate the mass loss rate and the terminal velocity of 
the wind to stellar parameters like mass, radius and luminosity. In our case, these 
parameters are the mass of the companion, the radius of the companion R^. and 
the luminosity of the companion L^. Some well-known stellar- wind models which 
we consider here a r e those due to Castor et a/., Kudritzi-Re i mers, and Vink al. 



fICastor et al. 1975 



Reimers 1975 



Kudritzki &: Reimers 1978 



Vink et al. 20001 ) . which 



span the range fro m the classic CAK model (jCastor et al. 19751 ) of the 1970s to the recent 



Vink et al. model fiVink et al. 20001 ). 



The rates of mass loss in these models are given by: 



Kudritzki-Reimers : M^ = 7kr 



LnRr 



Mr 



Castor et al: = 7cak 



Vink et al. M^ = 7vink 



MO-5 

^2.45^-0.5 



(14) 
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The constants 7kr, 7cak, and 7vink are given in the original works. We use the KR 
model in our main work here, and summarize the effects of varying the stellar-wind model 
in SecEXH 

5.1.1. Stellar models 

Results of numerical calculations of stellar evolution give stellar parameters like radius 
and luminosity as functions of the mass, metallicity and the age of the star. For using 
the above stellar-wind models in our calculations, we need these parameters for a given 
mass and metallicity of the wind-driving companion. In the spirit of our semi-analytic 
approach, we use the parameters given by Hurley, Pols and Tout (2000, henceforth HPT) in 
their comprehensive work on the construction of analytic approximations to the standard 
numerical results of stellar-evolutionary codes. The range of M^. shown in Fig|T] and 
discussed in Sec J4. 3. 11 leads to the conclusion that, except for a small region at the lower 
end of this range, the massive companions ignite He in the core while in the Hertzsprung 
gap (HG), while those in that small region do so at the top of the giant branch (GB). The 
companion's luminosity remains nearly constant in the HG, as discussed in HPT. Relevant 
fitting formulae for the luminosity and radius in the core Helium burning (CHeB) phase are 
given in Sec. 5 of HPT, particularly Sees. 5.1 and 5.3. We have used a metallicity z = 0.02 
in our main calculations, varying it later to test its effect on HMXB distribution (see 
Sec l6A3ll . 

The terminal velocity Vw which appears in Eqn{T3] is of the order of escape velocity 
from the surface of the companion and hence can be taken as a function of Mc, so that the 
accretion rate can be expressed in the form: 
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M^l^, (15) 

where the function f{Mc) ~ M^/V^ contains the physics of wind mass-loss. Equation 
gives the accretion rate as a function of Mc and a and so, along with EqnJTT] provides 
the transformation from (M^, a) to (M,Pi,), 



5.2. Transformation of parameters 

Inverse transformations for eqns. [15] and [11] need to be obtained to transform the 
PDF. We first carry out the transformation from (M^, a) to (M, Pf,). We note here that the 
PDF as a function of M and L will be identical up to a numerical factor since the two are 
linearly related to each other. The inverse transformation can become complicated because 
f{Mc) can be a very complicated function and hence difficult to invert analytically. Our 
goal would, of course, be to construct a procedure that works with a general /(M^), so that 
the formalism can be applied to any model without changing the procedure. We first write 
the accretion rate as 

M_io = 3.92^^^^. (16) 

The numerical factor is so chosen that 7kr = 10"'^; for other models the numerical 
constants can be appropriately adjusted in a similar way. With the aid of Eqns{TT] and [TS] 
we can eliminate a and write 

g{M,) = G{P,,M.,o), (17) 



where 
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(M, + Mns 



2 ' 

G{Pb,M.,o) = 2.78 X 10-^M!ioP,^ (18) 

This equation can be numerically solved to obtain as a function of {Pb,M-io), 
which is the first inverse transformation equation. After computing Mc numerically, a can 
be computed using 

a = h{M,)H{P,,M.io), (19) 

where 



H{Pk, M-io) = 3.3 X 10-2m_ioP^ (20) 

Using these inverse transformations, we calculate the Jacobian. The partial derivatives 
take the following forms in terms of the functions G{M_iq, Pf,), H{M_iq, Pf,), g{Mc) and 

dMc dG f dg 



dM^ dG f dg 



dPb dPb \dM^ 

da , dH dh f dq \ dG 

h — : VH ' * 



dM_io dM_io dM^ \dMj dM_io 

da , dH dh f dq \ dG , , 

h^Tf^ + H -— (-^] — (21) 



dPb dPb dM^ \dMj dPt, 

From eqs. f|T7|) . f|T9|) & fl2T|) . we can calculate the Jacobian as a function of Pb and M_io 



as: 
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J 



/(Me) 




3 



df 



2 



-1 



1.84 X 10"^ M^ioPb 



(22) 



5.3. 



Transformation of Distributions 



The distribution of HMXBs in Pb and M_io can now be obtained witli tlie aid of 



Eqs.dH]) and It is: 



fHMXBiPb, M-IO) = /P5JV (Me(P6, M-lo), «(n, M-lo)) | M_io)| , 



(23) 



and M and L are of course proportional to each other, as explained above. 



The resultant bivariate PDF fuMXE^L, Ph) is displayed as a surface in FigJTJ and the 
individual PDFs for the L- and Pfc-distributions are shown in FigsO and [S] respectively. For 
obtaining the individual distribution of each variable, the bivariate distribution is integrated 
over the other variable, as explained earlier, the range of integration extending over only 
the allowed range of the concerned variable. 

In Figjni we show the allowed region in the L — Pb plane, which shows how the allowed 
region in the Mc — a plane shown in Figj5] is mapped onto the plane of the new variables. 
The effect of the above transformation is a rotation of the allowed zone in {Pb, L) plane. 
The simple bounds on parameters of primordial binaries given earlier transformed into a 
nearly rectangular allowed region in the parameter space (Mc, a) for post-SN binaries (see 
Fig|5]). After transformation to the HMXB parameters, however, this allowed zone rotates 
into an inclined band, as shown in FigJHl The boundaries given in this figure serve as 
integration limits for the calculation of the individual PDFs from the bivariate distribution 
in terms of L and Pb- 



- 28 - 

6. The HMXB Distribution 

We now discuss the nature of the distributions of the luminosities and binary periods 
of HMXBs calculated above, and compare them with the current state of observational 
knowledge of these distributions. We emphasize that we are not attempting a detailed fit 
to the data at this stage, but focusing instead on a comparison between the general trends 
in calculated and observed distributions. The purpose of such a comparison is of course to 
determine if a more detailed computational scheme, e.g., a numerical population synthesis, 
is worthwhile in the future for a giving a more detailed account of the observations. A 
major virtue of an approach like ours is that it is capable, at least in principle, of assessing 
in a transparent way the relative importance of various components {e.g., the role of the 
primordial binary distributions vis-a-vis that of the evolutionary processes like the first 
mass transfer and the SN) in determining the final HMXB distribution. Such an assessment 
helps greatly in planning the strategy of subsequent computational studies. 



6.1. Orbital Period Distribution 

The PDF of HMXBs as a function of Pb is obtained immediately by integrating Jhaixb 
given in Eqn|23]over the accretion rate. The computed Pb distribution without SN-kicks 
(see fig. [HD is rather similar to the corresponding post-SN a distribution, i.e., a fiat Opik's 
law in the mid- region {Pb ~ 300 — 3 x 10^ hrs), with a gradual rise to this flatness at 
shorter periods and a slightly sharper fall-off from it at longer periods. The distribution 
with the inclusion of SN-kicks also resembles the post-SN a distribution with SN-kicks: the 
rise now is more gradual, reaching a maximum at Pb ~ 1000 — 1500 hrs. Instead of a fiat 
top, a slow fall is observed mid-P;, range, followed by a sharper fall-off at long periods, as 
before. The major and obvious difference between the two distributions is an overall shift 
towards wider orbits and longer periods when SN-kicks are included, entirely as expected 
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and as seen earlier in the post-SN a-distribution. We compare our theoretically obtained 
Pfe-distribution with observations in Sec j6.3[ 



6.2. X-ray luminosity function 

The XLF of HMXBs can be obtained by integrating Eqnl23] over Pf,, using the limits 
described previously and carrying out a simple linear transformation for M_io to L^q. The 
XLF is given by: 

hiLse) = / fHMXBiPb, M-io) dP, (24) 

where M_io = L^q/LIT. Figure |9] shows the numerically computed XLF. A broken 
power-law can be fitted to this computed XLF, with a cut-off at the neutron-star Eddington 
limit. Equation [23 gives the power-law exponents in the low and high luminosity regimes, 
obtained by least-squares fits to the computed XLF: 



dN -1.63 0.65 < L36 < 7.5 

^oc with r]= { (25) 

-0.98 10-3 < L36 < 10-2 



dL 



The XLF calculated without SN-kicks is also given for reference. It overlaps with the 
XLF with SN-kicks in the low-luminosity regime, and shows a slightly shallower power-law 
of exponent is rj = — L43 in the high-luminosity regime, the cross-over point being at 
Lcr ~ 8 X 10^^ erg/s. 

We now examine the effects of the various model parameters, e.g., stellar-wind model, 
IMF slope, metallicity, and so on on our theoretical XLF. 



- 30 - 



6.3. Comparison with observations 
6.3.1. The XLF 

We first compare the calculated L-distribution of Fig. [9] with the observed distribution. 
The major feature of the observed L-distribution is a power law with a differential slope 
dN/dL oc L~^'^ over a wide range of luminosities 3 x lO^^erg/s < L < 3 x 10^°erg/s. 
This is the so-called "universal" X-ray luminosity function (XLF) of HMXBs, obtained 
in the following way. Whereas the observed XLFs of various nearby early-type galaxies 
{e.g., the Milky Way, SMC, M82 and M83, the Antennae, NGC 4736, and so on) follow 
this trend, their normalizations are not the same. However, when the XLF of a given 
galaxy is normalized by the current star-formation rate (SFR) in that galaxy, XLFs for 



all th ese galaxies fa 



XLF ( 



Gilfanov 2004 



1 essentially on top of each other, yielding t 



le above "universal'' 



Gilfanov et al. 2004a 



Gilfanov et al. 2004b 



Grimm et al. 2002 



Grimm et al. 20031 ). This is due to the well-known fact that the strength of the HMXB 
population and their X-ray output is proportional to the current SFR, and is closely related 
to the discussion given in earlier on how the rapid HMXB evolution leads to a distribution 
of HMXBs which is, in effect, a "snapshot" of the galaxy taken at its current SFR. 

In comparing calculated and observed XLFs, we note first that, since we confine 
ourselves to only neutron-star HMXBs in this work, our calculated XLF applies only to 
luminosities not exceeding the Eddington luminosity for a ~ I.AMq neutron star, i.e., 
~ 2 X lO^^erg/s. The observed XLF extends upto almost 2 decades of luminosity above 
this, and the HMXBs at these higher luminosities are believed to be black-hole systems, 
with the possibility of both stellar-mass black holes and intermediate-mass black holes (at 
the highest luminosities) being present. We return to this question in more detail in Sec J7.3[ 

With this caveat in mind, we note that there is a remarkable agreement between the 
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calculated and observed XLFs in the luminosity range 3 x lO^^erg/s < L < lO^^erg/s, above 
which the calculated XLF cuts off at the Eddington limit for neutron stars, as expected, so 
that no comparison with observations is possible at higher L. Below L ^ 3 x lO^^erg/s, 
there was no data when the works referred to at the beginning of this subsection were 
published. However, more recent XMM — Newton observations of the Magellanic Clouds 



have extended the XLF be 



ow this lower limit, down to about L ^ 10^ erg/s or slightly 



lower, for SMC and LMC flShtvkovskiv fc Gilfanov 2005al : 



Shtvkovskiv fc Gilfanov 2005bh . 



The results indeed suggest a flattening of the XLF at low luminosities, and the best-fit 
differential slope of the observed XLF at these luminosities for SMC, which is given in 
the above reference as — 1.13lo''^3, ^^^^ almost consistent with our calculated slope 
at low luminosities, given in Figj9]and Eq. (!25|) . However, we must be cautious with our 
calculated XLF in the low-luminosity regime, since we have neglected (a) main-sequence 
companions undergoing atmospheric Roche-lobe overflow, and (b) Be-star companions in 
our calculations here, as explained above. 



6.3.2. The P^- distribution 

Observations of HMXB orbital periods are available in substantial numbers only for our 
galaxy and the M agellanic Clouds . Measurements of orbital periods have been performed 



for ~ 80 systems (iLiu et al. 2005 



Liu et al. 20061 ). Figur ^TO] shows the distribution of 



orbital periods constructed with this data. The distribution shows a rising part in the range 
10 — 300 hrs. With the (large) error bars, the distribution in the region Pb ~ 300 — 10000 hrs 
is consistent with either a uniform trend or a slow rise/decay with a peak around Pf, ~ 2000 
hrs. It can be seen that our theoretical P^-distribution is generally consistent with the 
observed one within the error bars, except at very long orbital periods, where the observed 
distribution shows a sharp cut-off beyond orbital periods ~ 1 year. Such an apparent cut-off 
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is expected for two observational reasons. First, it is very difficult to follow systems with 
such long orbital periods for sufficient times to establish reliable orbital periods. Second, 
such wide binary systems would generally have low luminosities, which would put them at 
the faint end of the XLF, and add further to the difficulties of a successful observation. 
Thus strong selection effects work against observation of HMXBs with long orbital periods. 
By contrast, the theoretical distribution extends upto the widest orbits that can survive 
the basic processes involved, particularly the SN explosion. The fact that these very wide 
binaries do not show up in the observed distribution does not imply that they do not exist, 
but simply that they are unobservable in practice. 

6.4. Parameter Study 

We now examine the effects of the various model parameters, e.g., stellar-wind model, 
IMF, metallicity, and so on on our theoretical XLF. 

6.4.1. Effects of stellar wind models 

We consider the two alternative models by CAK and Vink et ai, which were introduced 
in Sec l5.ll XLFs computed for these models are shown in FigJTTl It can be seen that the 
XLFs overlap in the low-luminosity regime for all models, which means that the XLF is 
independent of the exact form of /(M^) (see above) in this region. The XLF's differential 
slope is close to -1 in this regime. The position of the cross-over luminosity Lcr is different 
for the three models with similar power-law exponents above respective Lcr- The CAK 
model cuts off rather abruptly at luminosities considerably below the Eddington luminosity 
of neutron stars and so appears unable to account for the observations. Barring the 
difference in the Lcr, the Kudritzki-Reimers model and the Vink model give similar results. 
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in general agreement with observations. 

6.4.2. Effects of the IMF slope 

The IMF is typically given as a single power-law in the mass range that is relevant to 
the problem at hand. The power-law index of the IMF (a) is considered to be within the 
range ~ 2.0 — 2.7, a = 2.35 being the standard value given by Salpeter and widely used 
in this mass range. We study the effect of varying a on XLF in the range 2.0 - 3.0. The 
effect is shown in the FigJTTl The XLF at low luminosities below the kink shows hardly any 
variation. We therefore display the XLF only for luminosities above the kink. Even above 
the kink, the XLF slope varies by a small amount (from -1.47 to -1.66 between the lower 
and upper limits to a considered here), showing relative insensitivity of the XLF to the 
IMF slope. For an exact matching with the observed XLF, a slightly steeper IMF seems to 
be preferred, if other parameters remain unchanged. 

6.4.3. Other effects 

Among other parameters in the problem, we consider the stellar-wind velocity and 
the metallicity z. The terminal velocity of the wind is an important factor in determining 
the capture fraction. It is typically of the order of the escape velocity at the surface 
of the companion and generally thought to be ~ 10^ km/s. We varied around this 
canonical value to test the effect on the XLF. The changes in both the XLF slope and the 
high- luminosity cutoff were insignificant. 

Consider next the metallicity 2;, which affects the stellar parameters and the wind mass 
loss rate. We considered a large range of z from 10~^ to 0.02, and used the model of Vink 
et al. (2001) to study the effects of varying z on our final results. The results are shown 
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in FigJT2l There is no significant difference between the XLFs for z = 0.01 and z = 0.02: 
both have a slope ~ —1.6 in the luminosity range 2 x 10^^ — 2 x 10^^ erg s~^. The XLF for 
z = 0.001 also has nearly the same slope, but a somewhat different crossover luminosity. 
Finally, the XLF for z = 10"'^ is markedly different, but such low metallicities are not 
realistic for HMXBs, as they would be expected only in old stellar systems. 



7. Discussion 

In this work, we have developed a strightforward, first-principles scheme for 
understanding collective properties of HMXB populations, starting from standard, 
well-known collective properties of primordial binaries which are the progenitors of HMXBs, 
and following the transformations of the probability distributions through the evolutionary 
processes that lead to the formation of HMXBs. Our purpose, of course, was to assess if the 
standard picture of primordial binaries and the standard evolutionary scenario for HMXBs 
together can account for the observed collective properties of HMXBs in a basic, simple 
way. The fact that we find that such an account can indeed be given is most encouraging, 
and it constitutes, in our view, an essential step towards attempts at elaborate population 
synthesis schemes designed for understanding further details of these collective properties. 
Since our procedure is transparent and readily understandable at each step, it is easy in our 
scheme to follow the role of each ingredient in shaping the final HMXB distribution. In this 
section, we discuss various issues of principle and procedure which are relevant to this line 
of approach, and conclude with plans for the future. 
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7.1. Primordial-binary mass distribution 



We pointed out in Sec J2.1l that distributions of the masses Mp, Ms of the primary 
and secondary in the primordial binary can be described either in terms of these 
masses themselves, i.e., the pair {Mp,Ms), or alternatively in terms of the primary 



mass and the m ass ratio q 



M.jMp, i.e., the pair 



Mp,Q). Both have been done 



Jaschek fc Ferrer 1972 



Halbwachs 1983 



Trimble 1990 



in the literature ( Warner 1961 

^ I 

Kouwenhoven et al. 2007i : iKobulnicky fc Fryer 2007i ) , and we have chosen here the second 

description for its closer correspondence with essentially all recent work. The relation 
between these two approaches has been discussed thoroughly by Tout (1991). 



7.2. Be-star binaries 

Massive companions in HMXBs are of two types in general, viz., OB giants/supergiants 
and Be stars, the latter being characterized by (a) strong, broad emission lines that supply 
evidence for rapid stellar rotation, and (b) somewhat lower masses and wider orbits. Since 
we have confined our detailed calculations in this work to the former, as stated in SecJS], we 
now discuss the expected role of Be-star binaries in HMXB XLF. 

It is believed that Be stars are often surrounded by an outflowing disk of matter 
expelled by centrifugal forces from the fast-roating equatorial regions of the star, in addition 
to the usual stellar wind emitted from all over the stellar surface. Accretion by the neutron 
star from this outflow material generally follows the basic description from stellar winds 
given in Sec J5.1[ with appropriate values of for the fast wind and the slowly outflowing 
disk, with one major caveat. Since the matter in the rapidly-rotating outflow disk has much 
angular momentum, it may form an accretion disk around the neutron star, which would 
then drain on the neutron star a slow, viscous timescale. This complicates the description 
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considerably, as the outflo w disk is generally expected to be tilted with respect to the 



orbital plane ( IGhosh 1995[ ). so that the orbiting neutron star would "crash" through this 
outflow disk (twice per orbit in general), possibly acquire an accretion disk, and accrete it 
slowly over the rest of the orbit. This would naturally lead to outbursts of X-ray emission, 
which are indeed observed in Be-star HMXBs. 

The question for our purposes here is: how does all this affect the XLF of HMXBs? 
First consider the observational situation. The point to note here is that Be-star HMXBs 
are basically transient systems with low duty cycles, so that a long-term average of the 
luminosity of a Be-system is much lower than that of an OB supergiant-system. Thus, in an 
XLF constructed from a long-term monitoring of the X-ray sky with an all-sky monitor, one 
would expect the high-L parts to be dominated by OB-systems, while the low-L parts may 
have considerable contributions from Be-systems. Indeed, the observed XLFs cited earlier in 
this paper have been constructed from recent work with X-ray observatories in the following 
way. The XLF for the Milky Way has be constructed from ~ 4 years of observation with the 
all-sky monitor on RXTE. Thus, the above argument certainly applies to this case. Indeed, 
since the typical long-term average expected from Be-systems would be at or below the 
lower end of the luminosity-range over which the XLF is actually reported (and this applies 
to basically all observed XLFs except those obtained from XMM — Newton observations 
of SMC and LMC; see Sec J6.3.1|) . we would expect little contribution to the reported XLF 
from Be-star systems. 

However, note that the XLF for other nearby galaxies reported in these references 



have been constructed essentially 



Tom one "sna pshot" {i.e., a single exposure) taken by 



Chandra and XMM — Newton (IGilfanov 20041 ). Given this, it is remarkable that the 



XLFs of all these galaxies (suitably normalized by their SFRs) are essentially coincident 
with one another over the same range of L, as explained earlier. Considering the fact that. 
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because of their low duty cycles, only a small fraction of the Be-systems would be present 
in these "snapshot" XLFs, we would still expect some contribution from them in a range 
of L typical of the outbursts of Be-systems. The fact that these XLFs appear very similar 
to the above long-term average XLF of the Milky Way over the canonical luminosity range 
3 X lO^^ergs^^ < L < 3 x 10'*°ergs~^ is therefore most noteworthy, and may imply that, 
for reasons which are not clear at present, Be-systems may not have made a substantial 
contribution to these observed XLFs. 

Now consider the inclusion of Be-systems in calculational schemes like ours. To the 
extent that the wind- accretion formalism can be applied to accretion from both the fast, 
low-d ensity wind 



Tom the stellar surface and the slow, high-density outflow in the equatorial 
disk ( iGhosh 19951 ). these systems are already in our scheme, at least in principle. However, 
while a quantitative description of the periodic acquisition an d drainage of accretion disks 
described above has been done for individual binary systems (jPravdo &: Ghosh 200ll ). its 
inclusion in a study of collective properties of HMXB populations is more complex and 
outside the scope of this work. 



7.3. Black-hole HMXBs and ULXs 

As mentioned in Sec J6. 3.1] the observed "universal" XLF of HMXBs extends to about 
two decades of luminosity above the Eddington luminosity for canonical IAMq neutron 
stars. While we focus on neutron-star HMXBs in this study, it is interesting to consider 
this brightest end of the XLF briefly. If the X-ray sources here are accretion-powered, they 
can only contain accreting black holes (some sources upto luminosities ~ 10^^ erg s~^ can 
of course be close juxtapositions of several neutron-star sources, as has sometimes been 
suggested), either stellar-mass ones {Mbh ~ IOMq), or even the so-called intermediate-mass 
black holes (IMBH), with masses Mbh ~ (10^ — 1O^)M0, the X-ray sources corresponding 
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to the latter objects being often called ultra-luminous X-ray sources (ULX). 

A truly remarkable feature of the universal HMXB XLF is that a single, smooth power 
law gives an excellent account of X-ray binaries containing neutron stars, stellar-mass black 
holes, and IMBHs over the entire l uminosity range 3 x lO^^ergs^^ < L < 3 x 10^°ergs~^, 



i.e., above the "kink" in the XLF (IGilfanov 20041 ) . While it is not difficult to extend 



the formation and evolution scenario outlined in the earlier sections to include higher 
primary masses which would produce stellar-mass black holes, and then imagine that the 
other systematics would go through in such a way as to extend the power-law XLF, this 
argument does not automatically include ULXs and IMBHs, whose formation scenario has 
to be different and more exotic, e.g., black- hole mergers in dense stellar clusters. And 
yet, as Gilfanov (2004) has noted, these "rare" and "exotic" objects appear to form a 
smooth extension of the "ordinary" HMXB population. Whether this is really true or not 
can possibly be probed with future observations in a way suggested by this author: if we 
adopt the alternative hypothesis that the apparent cutoff in the currently observed XLF at 
L ~ 3 X lO^^ergs"^ really corresponds to the maximum possible luminosity of what we may 
call "ordinary", stellar-mass black holes referred to above, then these "exotic" IMBHs may 
show up beyond this cutoff if we can observe regions with extremely high star-formation 
rates, since the merger scenario for the formation of IMBHs implies that they would occur 
in very dense regions with very high star-formation rates. However, since such IMBHs 
would necessarily be rarer than stellar-mass black holes, we should expect a "step down" at 
the presently observed cutoff, beyond which the XLF would continue at a lower strength. 
This is a fascinating possibility. 
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7.4. Shape of the XLF 

Understanding the XLF shape is an important step towards extracting information 
about various processes which determine the coUective properties of HMXBs. Two main 
features of the XLF calculated using our scheme are (1) A kink at Lcr ~ 8 x 10'^^ erg/s, and 
(2) a power-law behavior with the differential slope of -1.6 above Lcr- We will discuss the 
possible origin of these two features in this section. 

7.4.1. TheXLFkmk 

The origin of the kink in the HMXB XLF can be understood in terms of allowed zones 
in the parameter space. This is demonstrated in a clear way by overplotting contours of 
constant L on the allowed zone in the (M^, a) plane, as shown in FigjSl Let us first examine 
a simple "toy" model for the XLF. In this model, we approximate the post-SN PDF as a 
function of Mc and a as fpsn oc M'^'/a (as discussed earlier, this "toy" is not too bad an 
approximation over much of the parameter range in case of no SN-kicks). We now note that 
L can be written in a schematic way as L oc f{Mc)/a?. Next we transform the approximate 
form of the post-SN PDF from (Mc, a) to (Mc, L). A bit of simple algebra shows that such 
a transformation gives a PDF of the form fu oc M^'^/L. Integrating this PDF over Mc 
would yield dN/dL oc , if the limits of integration were independent of L. Figure O 
shows that this is indeed the case for lO^^erg/s < L < L^r, since the Mc-limits are clearly 
seen to be almost independent of L when L lies in this range. 

By contrast, as we move to higher L- values, the contours become shorter by cutting 
off low- Mc regions. In other words, the lower limit of integration becomes a function of L. 
Shortening of the contour length results in a power-law XLF with the slope steeper than 
-1. Eventually, near the neutron-star Eddington luminosity, the contour passes out of the 
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allowed region altogether, and the XLF is cut off. 

At very low luminosities we expect a similar mechanism to cause a turnover in the XLF 
(which would cut it off eventually at extremely low luminosities), as is clear from FigJSl 
since higher- Mc regions are progressively cut off as L decreases in this range. However, 
this effect is not expected to be very important for two reasons. First, cutting of the 
highest-mass regime in is not a severe problem, since the Mc-distribution drops steeply 
in that regime anyway. Therefore, one needs to go to very low luminosities (< 10^^ erg/s) 
to observe this effect. Second, with current observational sensitivity, we can track the XLF 
behavior below only for a few nearby galaxies, and the very low luminosities indicated 
above are not even approached. It appears, therefore, that this XLF turnover is unlikely to 
be amenable to observation in the near future, and accordingly we do not consider it any 
further. 



7.4.2. The XLF slope 

The XLF of HMXBs obtained observationally does not reach luminosities below L^r 
for most of the galaxies (see Sec J6.3p . Therefore, only the power-law regime above Lcr can 
be compared with observations in most cases. Physical origins of this power-law index are 
an important aspect of our understanding of the collective properties of HMXBs. Note 
first that our calculation of the XLF is a numerical one, and the power-law result given in 
Sec. 16.21 is only an analytic approximation to it. However, simple qualitative arguments 
may serve to illustrate the basic physics underlying such calculations, and we consider such 
arguments below. 

In a pioneering argument of this type, Postnov (2003; also see Postnov & Kuranov 
2005) proceeded as follows. Expressing the XLF as: 
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dN _ dN dM^ 

'dL~dMc~dL' ^ ' 

one can use a suitable estimate of the mass- function dN/dMc of the companion for 
these arguments. The estimation of is more involved, and the original Postnov method 
was to estimate as follows, and use its reciprocal. Using L oc M oc M^/ (a^V^) (see 
SecjS]), this author obtained a power-law form L oc , with the aid of a simple stellar- wind 
model oc Lc/Vu, oc Lc^/Rc/Mc, into which the following scalings for massive stars were 
inserted: Lc oc Mj}, Rc oc M°'^. This gave /3 ~ 2.5. Finally, assuming a power-law form 
for the companion mass distribution, i.e., dN/dM^ oc M~", straightforward algebra with 
Eg. (1261) leads to an XLF slope of dN/dL oc 1-^"+^-^'^/'^ . Postnov (2003) assumed a = 2.35 
{i.e., the exponent of the Salpeter IMF) for the Mc-distribution, which led to an XLF slope 
of ^ -1.54. 

In revisiting the above argument, we note first that Eg. (1261) is incomplete, since N is a 
function of both Mc and a, whose form we have calculated explicitly in the earlier sections. 
Thus, the complete eguation is 



dN _ dN dMc dN da 

11 ~ dWc'dL ^ ~dadL ' ^ ' 

Next, we note that the scalings of Lc and Rc with Mc used in the Postnov (2003) work 
apply to massive main-sequence stars, but not to the evolved massive stars of interest here. 
Also, the wind mass-loss prescription used in that work is similar to that in the CAK model 
of the 1970s, which our calculations have shown to be inadeguate. The rough scalings 
for Lc and Rc for evolved massive companions of interest here can be obtained from the 
appropriate formulae in Sec. 5 of HPT, and are Lc oc M^-^, Rc oc M^-®. 

With these scalings, and our Mc-distribution shown in Fig{T], which we can approximate 
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with the power-law of exponent a = 2.4 that apphes to its principal part, we can 
evaluate the first term on the right-hand side of Eq. fl27l) . For the Kudritzki-Reimers 
model, it can be easily shown that /3 ~ 5.4, so that this first term gives an XLF slope 
— {a + /3 — ~ —1.27. The other stellar model gives an essentially identical final result. 
The rest of the contribution comes from the second term on the right-hand side of Eq. (!27l) . 
which can be easily calculated and which leads to our overall XLF slope ~ —1.63. 

We have gone through this argument in detail because it addresses an interesting 
observation made by Gilfanov (2004) that, in HMXB systems powered by stellar-wind 
accretion, the distribution of M and therefore L should be governed by the properties of 
the massive companion, in particular the distribution of Mc and Lc. Our work here shows 
that this is largely, but not completely, true. The properies of the massive companion 
are contained in the first term on the right-hand side of Eq. (|271) . and keeping only this 
term, as Postnov (2003) did, amounts to neglecting altogether the binary orbital properties 
which are contained in the second term on the right-hand side of Eq. fl271) . and which also 
influence the XLF. The estimates summarized above give a measure of the relative sizes 
of the effects of the companion and the orbit, and demonstrate that, while the latter are 
certainly smaller, they are by no means negligible. 



7.5. Conclusions and outlook 

In this work, we have described a method for obtaining the distributions of some of the 
essential collective properties of HMXB populations in the stellar fields of normal/starburst 
galaxies, wherein we start from accepted distributions of primordial binaries which are 
progenitors of such HMXBs, and follow the transformation of these distributions with the 
aid of a Jacobian formalism as the primordial binary population evolves into the HMXB 
population. Our method, which is semi-analytic, traces in a transparent way the effects of 
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various processes in the course of this evolution, and so assesses with ease which physical 
processes dominate in determining which distribution. For example, the distribution of 
the properties of the massive companions seems to have the dominant effect on the XLF, 
although the distribution of the orbital parameters does have a significant effect, as we 
demonstrated in Sec J7.4.2l But the distribution of HMXB orbital periods appears to be 
strongly influenced by both the primordial orbital distribution and the SN-kick properties. 

The agreement between our calculated XLF and binary-period distribution and the 
observed HMXB distributions is most encouraging, and it justifies a future Monte Carlo 
population synthesis scheme for a more detailed undrstanding of how HMXB populations 
are built in the stellar fields of normal galaxies. However, we must first extend our present 
method to the more complex problem of following the formation and evolution of LMXB 
populations from their corresponding primordial binaries, as mentioned in SecJH This will 
occupy us in the next papers in this series. 

It is a pleasure to thank M. Gilfanov, E. P. J. van den Heuvel, V. Kalogera, L. Stella, 
and R. A. Sunyaev for stimulating discussions, and an anonymous referee for valuable 
comments which greatly improved the paper. 



A. Averaging supernova kicks and post-SN transformation equations 



A.l. Calculation of SN-kick average 



In this appendix we describe the method of averaging the effects of SN-kicks for a 
Maxwellian kick-distribution. It has been widely assumed t hat the distributio n of kicks will 



be isotropic and recent observations by have supported it (iHobbs et al. 20031 ). Therefore 



the angle dependence in Eql5]is averaged in a straightforward way, which leads to Eql6l 
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Two points are to be noted here. First, the only surviving kick-term in the averaged 
transformation equations is vf, which is to be obtained by averaging v"^ over the Maxwelhan 
distribution. Second, it must be remembered that the process of averaging gives the 
average values of post-SN quantities. These therefore should be compared with the average 
behaviour of a number of sample observational sets. 

A further, most crucial, point to be noted at this stage is that in our study here (and 
in all similar studies), we are interested in only those post-SN systems which remain bound 
as XRBs, so that they can eventually produce HMXBs whose distribution we are interested 
in. Accordingly, we must exclude at this point all systems which become unbound in the 
SN. The way we do so is as follows. We note that, in the f ^-averaging process, we would 
make an error if we carried out the integrations over the Maxwellians in Eq|7]upto infinitely 
large values of f , since this would include all the unbound systems that are to be excluded. 
Rather, we must truncate the integration at a suitable upper limit Vup which corresponds 
to the point at which the post-SN system becomes just unbound. This point is readily 
obtained from the first of Eqnsj6]by setting a/ — > oo, and Vup given by: 



We thus work with a Maxwellian which is truncated from above at v^p (and suitably 
normahsed), so that vl is given by: 



dispersion a. This simplifies the expression for f^, which can be written as = a'^h{f), 




(Al) 



v^exp{-vy2a^)dv 
JJ'"'' v'^exp{—v^ /2a^)dv 



(A2) 



Vup is well-defined for each set of values of pre-SN parameters. For algebraic 
convenience, we can express it as Vup = fy/2a for a given Maxwellian distribution with 
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where h{f) is given by: 



h{f) 



3- 



hi{f)-h2{f) 



h2{f) 



hiif) 



(A3) 



FigJT3] shows h{f) as a function of /. Two regimes are clearly demarcated, with a small 
transition region between the two. The low-/ region can be described by a power-law given 
by, h{f) = 1.2/^, whereas in the high-/ region, h{f) = 3, which is its asymptotic value. 
The crossover point of the two regimes is at fc ~ V2.5. This clear division in two regions 
give us a straightforward relation for f^. For / < f^, vl is independent of a, and is given 
by vl = O.Gvlp. On the other hand, for f > fc vl is independent of v^p, and is given by 
vl = 3cr^. This behavior of h{f) can be understood as follows. For / < fc, Vup is small and 
hence only the initial rising part of the Maxwellian is relevant. Thus the distribution is 
essentially given by f{v) oc v^, which is independent of a. On the contrary, for / > fc, Vup 
is large and almost the entire Maxwellian profile is included, excluding only a small tail. 
Detailed calculations given above show that the transition region between these two regimes 
is small, so that, as a first approximation, we can use the following simple prescription: 



The above prescription can be substituted appropriately in eqn. [6] to obtain the 
post-SN parameter transformation or the inverse transformation, which is required for the 
Jacobian formalism. We describe in the next section the relevant transformation equations 
and the Jacobian for the two cases. 




f<fc 



f>fc 




(A4) 
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A. 2. Post-SN parameter transformation 



For f>f^ 



In this region, is independent of the parameters of the binary system. We can 
simplify Eqnj6]by using following ratios. Let rj = M*/aj, Vf = Mj/af, and = cr^. Also, 
we define orbit size change factor as s = af/ai. Equation [HI can then be used to write the 
inverse transformations: 



These equations are indirect transformations which give and M* on substituting the 
definitions of s and in above equations. Since the companion mass is unchanged in the 
SN explosion, one can readily calculate Mp ^ from M* for known M^. The Jacobian of this 
transformation can be calculated as follows. Note first that the transformation is only in 
two variables i.e. {af,e) — )■ (aj,Mpc), since the third parameter i.e. is unchanged. For 
the same reason, one can also write dMp^c/d{a,e) = dM^ /d{a,e). All the necessary partial 
derivatives can be written in terms of the derivatives of s and as follows: 



s 



1 + a/c^ (1 + Tk/rf) - rk/rf 
1 -e2 

r/(2s - 1) - 3rfc 



(A5) 



1 ttf ds 



daf 
dai 
de 



s s"^ daf 
ttf ds 
de 



dMl 



d{af,e) 



^d{af,e) ^d{af,e) 
dMp^c dai dMp^c dai 
de daf daf de 



sn 



(A6) 
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For f < 

vl in this region is independent of a and written only in terms of Vup, which is given by 
eqn. lAll The inverse transformations in this case are given by: 



1 + V(3e2-l)/2 



e2 



M* = _£ [5s (1 - e^) - 4] (A7) 

where s is the orbital-size factor defined similarly as aj/ai. The Jacobian in this case 
is given by an even simpler relation. One can easily see that M* is independent of a/ 
and depends only on e. Therefore dM^/daf = 0. The Jacobian in that case is given by 
Jsn = {dMll de){dai / da f) . Straightforward algebra shows that 



5 M}e 

■^^^ = Yss{l-e^) + l ^^^^ 

These two inverse transformations can be applied in the relevant regions to transform 
the PDF with the aid of the Jacobian formalism. 
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Fig. 1. — Distribution of companion masses in post-SN systems. 
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Fig. 3. — Distribution of eccentricities of post-SN systems. 
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Fig. 6. — Allowed zone in L—P(, space shown covered by crosses. Upper and lower boundaries 
of zone as indicated. 
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Fig. 8. — Calculated distribution of HMXB orbital periods. 
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Fig. 9.— Calculated XLF of HMXBs. 
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Fig. 10. — Observed distribution of HMXB orbital periods. Siiown is the data for ~ 80 
systems in the Milky way, SMC and LMC. Prom the catalogues by Liu et al.: see text. 
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Fig. 12. — Dependence of XLF on model parameters: metallicity as labeled. 
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Fig. 13. — Dependence of h(f) on f 



